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Abstract 

We extend the volume of fluid method for the computation of two-phase 
flow to a higher order accurate method in two dimensions. The interface 
reconstruction by the PLIC method is thereby replaced by a periodic in- 
terface reconstruction. The advection step is reformulated and extended to 
higher order in order to account for the present interface representation. This 
periodic interface reconstruction describes the interface in terms of higher or- 
der periodic B-splines. Numerical tests verify that the theoretical order of 
convergence is indeed exhibited by the present method. 
Key words: B-Splines, Two-phase Flow, VOF, High Order Accuracy 

1. Introduction 

Two-phase flow can be found in many industrial applications. A popular 
method for the computation of two-phase flow is the volume of fluid method 



(VOF) [2fl, |2|. 



The volume fraction, the central object of the volume of fluid method, de- 
notes the ratio of the volume (area in 2D) occupied by one phase in a cell of 
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the computational domain to the cell volume (cell area in 2D). The volume of 
fluid method can be subdivided into two steps: the interface reconstruction 
step and the advection step. The interface reconstruction step computes the 
interface position at time t using the volume fraction field at time t. The 
advection step advects the volume fraction field from time t to time t + At 
using the reconstructed interface at time t. The volume of fluid method has 

n n 

its origins in the works of [7] and [14j . Substantial improvement of the in- 
terface reconstruction has been achieved with the piecewise linear interface 
computation (PLIC) method by Youngs 26] in 1982. However, the resulting 
interface is approximated by piecewise straight lines which makes it necessary 
to estimate the curvature by additional approximation schemes, for example 



the height function method 



10 



9. 



11] . In order to obtain a more smooth 



interface, Price et al. in 1998 [17J derived a method replacing the straight 
lines by parabolas. Due to the fact that a numerical minimization has to be 



30 



a, the 



12| used 



performed in each cell to find all the coefficients of the interface para 
method enjoys less popularity. More recently, in 2004, Lopez et al. 
a parametric cubic spline interpolation through the midpoints of the PLIC 
interface lines and obtained a smoother description of the interface. How- 
ever, although cubic splines are known to interpolate a function with fourth 
order accuracy, their method inherits the second order accuracy of the PLIC 
method for the test cases presented in [12] since it is based on the same ap- 
proach. A further development of this interface reconstruction using splines 
to improve the interface obtained by the piecewise lines of the PLIC method 



has been presented in 



using quadratic splines. Both approaches are, how- 



ever, based on the PLIC method for reconstruction and advection and share 
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therefore also the drawbacks of the PLIC method. Anothe r appro ach replac 
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24j . where the 



ing this time the PLIC method has been presented in 
interface is divided into segments and each segment is reconstructed globally. 
This allowed for a more accurate description of the interface. A drawback 
of this method is however the need to choose a division of the interface into 
segments. 



In the present discussion we shall modify the approach presented in 
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22 



24| by deriving a periodic description of the interface separating two immis- 



cib 
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e li quid s in two dimensions. The present method is, as the method in 
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24| , a global method opposed to the PLIC method which uses only 
local information. The interface in the present discussion is represented in- 
directly by two functions depending on a periodic parameter. The actual 
position or other quantities, such as the normal or the local curvature at the 
interface, are then derived from these two functions. The advection step is 
adapted to the present interface representation. 

The present discussion is organized as follows: The present interface repre- 
sentation is derived in the next section, section [2j Periodic B-splines are used 
to approximate the interface, cf. section |3j In section HI the advection step 
is presented. The numerical verification is done in section Finally, the 
present discussion is concluded in section [6j 

2. A periodic representation of the interface 

In the present discussion we treat the case of a two dimensional drop 
of blue fluid enclosed in red fluid, cf. figure [TJ The red fluid occupies the 
domain Q re d, whereas the blue fluid occupies the domain Quue- These two 
domains are separated by a common boundary, the interface /. The central 



problem of the volume of fluid method is to compute the temporal evolution 
of the interface / when subjecting the fluids to a velocity field u: 



u y (x,y,t) 

Since we are dealing with incompressible fluids, the velocity is solenoidal. In 
the present discussion we assume that the interface / can be described by 
a periodic line I. We exclude topological changes in the present discussion. 
In addition any third phase should not be present in order to avoid contact 
points. We also assume the line I to be sufficiently regular. As for polygons, 
cf. [lj], the area of a domain VL b i ue , enclosed by a line I, can be computed by 
means of a function F: 

F{x) = l -x, (2) 

where x is the position vector of a point. The divergence of ((2]) is unity, as 
can be verified straightforwardly. Having now a periodic parametrization of 
the line /: 

I : x{s) = f X ^ I , s G [0,2tt) (3) 

V y{*) J 

the area V of fi Wue can be expressed by: 

^Mue ^blue 

(4) 

where the periodicity of the line has, without loss of generality, been chosen 
to be In. It is, in addition, implied that the tangential on I points in counter 
clockwise direction. We now define two functions a(s), resp. f3(s) by: 

s 

a(s) '= \ J x ( s ')y\ s ') ~ x'{s)y{s') ds', (5) 

so 



0{s) := x{s)y{s). (6) 

The derivatives of ct(s), resp. /3(s) are then given by: 

a'(s) = ±(x(s)y'(s)-x'(s)y(s)) (7) 
f3'(s) = x'(s)y(s)+x(s)y'(s). (8) 

Since the position x(s) of a point on the interface is periodic in s, we con- 
clude that the derivatives a', resp. 0' are also periodic in s. The position 
(x(s),y(s)) of a point on the interface on the other hand can then be recov- 
ered by the following expressions: 

x 



\0' ~ a' 



x 



(9) 



Integrating equations ([9]) and ffTUj) with respect to s gives us then the final 
result: 



x[s 



= xoexp / a(s')ds', (11) 

J SO 

y(s) = ?/o exp / b(s')ds', (12) 

J S() 

where (xo, yo) is the position of the interface for s = Sq. In order for equations 
d!]) and (|T0|) to be well defined we have to choose a coordinate system having 
the region Q re d in the positive quadrant sufficiently far from the origin. The 
area V included by the interface is then given by: 

V = a(2n). (13) 

The strategy of the present method is to represent the interface of the drop 
by the functions a, resp. and to obtain the position of the interface by 



formulae [TTJ resp. [T2J A normal n on the interface is then given by: 




Figure 1: The two dimensional domain il with the red fluid occupying the region fl re d 
and the blue fluid occupying the region fltiue- The common boundary of fired and Qbiue 
separates both fluids and is called the interface. At each point of the interface a normal n 
can be defined. We define the normal to point into the red domain. In this case where a 
blue drop is enclosed by red fluid this implies that we transverse the interface of the drop 
in counter clockwise direction. 



3. Interpolation by periodic B-Splines 



As mentioned above, instead of representing the interface position directly 
by B-splines, as for instance done in 25] in the framework of front-tracking 



methods or in 



12| for the volume of fluid method, we represent the interface 



by the two functions a, resp. (3, equations ([5]) resp. ([6]), defined on the inter- 
val [0, 2tt]. We use a uniform discretization of the interval [0,27r], meaning 
that we choose N + 1 knots Sj G [0, 2ir]: 



2m 



0. 



(16) 



dividing the interval [0, 2n] into iV sections of equal length. Periodic B- 
splines can actually handle more flexible discretizations, which could be used 
to distribute points to regions of interest. However, in the present discussion 
we restrict us to the uniform case. Having now the knots s$, i — 0, . . . , N, 
the periodic basis spline Bf of order P is obtained recursively by, see for 
instance 
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151: 



Bf{s) = {^^)Br i (s)+( Si+P+1 S )BZ?{s) (17) 

\Si+P — Si J \Si + p + i — Sj + i J 

1, Si < X < S i+ i, 

B°(s) = " (18) 

I 0, otherwise. 

The basis spline Bf has finite support [si, s i+ p + i). Therefore representing 
the function f(s) to interpolate as a linear combination R(s) of the periodic 
basis splines Bf: 

N-P-l 

i=-p 

leads to a cyclic banddiagonal system with bandwidth P for the unknown 
coefficients Cj. In the present discussion we will only use odd order B-sp 
The function values fi for interpolation are then taken at the knots Sj 

fi = f(si). (20) 

7or the resulting interpolation we have the following bound, see for instance 
1^ | . If / G C p+1 ([0, 27r]) and if / and the interpolating spline R of order P 
are periodic on [0, 27r] , the following bound holds: 

\\f-R\\ L *<C 2 (P)h p+1 \\f\\ HP+ i, (21) 

where h = maxi^Ar.! (sj — Sj_i) and the constant C is given by: 



mes. 
3: 



C(P) = 2-t( p+1 )- x (^±1 + 2^ 



(22) 



This implies that if choosing B-splines of order P the interpolation will have 
an order of accuracy P + 1 with respect to the grid spacing. The cyclic 
banddiagonal system resulting from ffl9l) can be solved efficiently by means 



of a banddiagonal solver in combination with the Woodbury formula 



16| 



The function j3, equation fl6]), is periodic and can thus directly be inter- 
polated by periodic B-splines. However, the function a is not periodic but 
takes different values at the right and left boundary of the interval [0, 2ir]: 

a(0) = a(2vr) = V, (23) 

where V is the area of the drop. The derivative a' is, however, periodic. 
In order to use periodic B-splines to interpolate the function a, we define a 
periodic function a* by: 

a*(s) = a(s) - ¥-s, (24) 

Z7T 

which is then interpolated using periodic B-splines. Once we have an ap- 
proximation to the functions a, resp. (3, equations (jSJ) resp. (JS}, we evaluate 
the integrals 

s 

A(s ,s) := J a(s')ds (25) 

so 

s 

B(s ,s) := J b(s')ds\ (26) 

SO 

from equations (jUJ), resp. ffTUl) by Gaussian quadrature [3], in order to 
compute the position by means of equations (jlip . resp. ffl2l . This will, 
however, introduce additional numerical error. A consequence of this is that 
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the integrals in equations (I2"5"]) . resp. (12"6"]) . might numerically not evaluate to 
zero for So = and s = 2tt. Therefore we determine first the total quadrature 
error e a , resp. e& by 

2V-1 

e a = Si,8 i+1 ), (27) 

i=0 
N-l 

tb = ^2q(B,Si,s i+1 ), (28) 

i=0 

(29) 

where the symbol q(A, Sj, s^+i) means taking the numerical quadrature of the 
integral A, equation (125]) . from Sj to Sj+i. Since B-splines are discontinuous 
in the P th derivative across the knots Sj, we perform a Gaussian quadrature 
on each subinterval [sj,Sj + i]. In order to assure that the position, given by 
equations (ITT]), resp. (]T2]) is itself a periodic function of s we replace the 
argument a, resp. 6 of the integrals in equations (125]) . resp. ( 126]) by a* and 
h* defined the following way: 

a*(s) = a(s)-^, (30) 

271 

b*( S ) = b{8)-£. (31) 

(32) 

4. Advection step 

Once we are given an interface I(t) at time t represented by the two func- 
tions a t and (3 t , equations ([5]), resp. (EJ), we need to formulate an advection 
scheme which allows to compute the interface I(t + At) at time t + At. Be- 
fore going over to the actual derivation we introduce the notions of flux and 
fluxing regions. 

9 



In the present discussion we assume that the flux M P0tPl through a section 
with end points po = (xo,yo), resp. p\ = (xx,yx) can be computed for arbi- 
trary points po and p\. The flux M POtPl is given by 

t+At t+At 

J Q P o,pi(t') dt' = J ip(x 1 ,y 1 ,t') - ip(x ,y ,t')dt', (33) 
t t 
where Q is the volume flow and ip the stream function. In the present dis- 
cussion the flux M po , pi will be computed analytically, since we are given the 
analytical stream function if) for the benchmark tests in section [51 The flux 
M po )P1 has a geometrical interpretation, cf. figure [2j It can be seen as the 
signed area of the region of points passing through a line with end points po, 
resp. p\ from t to t + At, the fluxing region. This region is bounded by the 
line from po to p±, its image at t, when tracing the line back from t + At to t 
and the trajectories r , resp. T\ of the points p , resp. p\. In order to trace 
a point p = from t to t + At we have to solve the following differential 
equation: 

^ = «(£,<), (34) 

with initial condition xq, where u = (—d y ip, d x ip) T is the velocity field, which 
is for the present benchmark tests, cf. section given analytically. We 
solve equation ( 134)) by the classical four stage Runge-Kutta method. For 
the present advection scheme it is necessary to approximate the trajectory r 
from time t to t + At of a point p = Xq. This is done by means of Lagrange 
oolynomials on the Gauss Labatto Legendre (GLL) nodes, cf. for instance 
18j . on the interval [t, t + At]. If n > 2 is the number of nodes chosen, the 
interval [t, t + At] is divided by the n GLL nodes into n — 1 sections 5ti, 
i = 1, . . . , n — 1. Solving equation f )34|) successively for each point in time 
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tj = t + YH=i 8U, j = 0, . . . ,n — 1 will give us a set of interpolation points: 

(ti.xfo) = (x(t,-),y(t,-)) T ), (35) 
where x(t ) = ^o- The trajectory r can then be approximated by: 

n-l 

T : x{t) = ^2x(t j )L j (t), (36) 

j=0 

where Lj(t) is the j Lagrange polynomial. 

The advection step is sketched schematically in figure [31 Having the interface 
I(t) at time t represented by the functions a t , resp. f3 t , equations ([5]), resp. 
(Q, we choose a sequence of points on the interface I(t), as depicted in figure 
[31 Several criteria might be possible according to which the points might be 
chosen 4|. However, in the present discussion we take the points pj at the 
parameter values Sj, the nodes chosen for the discretization, equation ffl6|) . 
The points are computed by equations ffTTl) . resp. f[T2|) : 

Pj = xj = x T (sj) = (xis^^isj)) , j = 0,...,N-l. (37) 

It is known that during simulation the points on the interface can cluster in 
regions of the interface if always the same points on the interface are traced 
forward [3j. However, in the present discussion we will not treat this problem 
but instead focus on the general method itself. 

Now that we have chosen a sequence of points Pj{t) on the interface at time 
t, we trace Pj{t) forward in time from t to t + At by solving (I34p . as depicted 
in figure [31 giving us the point pj(t + At). In the following the point Pj(t) 
at time t will be written with a tilde pj to indicate that this quantity is at 
time t. Its image pj(t + At) at time t + At will be written without tilde pj. 
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Focusing now on two consecutive points pj and Pj+i at time t + A, we know 
that the area bounded by the interface I(t) at time t, by the trajectories tj 
and Tj+i of the points pj, resp. Pj+i, and the interface I(t + At) at time t + At 
must equal M PjPj+1 , as depicted in figure El The flux M PjiPj+1 can thus be 
written as a sum of integrals along the bounding lines of the fluxing region: 

M Pj ,P 3 +i = Sitj + S Tj - S Tj+1 + S It+At p (38) 
where S/ t . is the integral along the interface I(t) from Sj to Sj + i, 

Sj + l 

S It . = J F{s')-n{s')ds' (39) 

= at t (s j+1 ) - a t (sj), (40) 

since we are representing the interface I(t) by means of at and (3 t , cf. equation 
(E]), resp. The integral S T on a trajectory r is given by 

t+At t+At 

/-4 7% ¥b n 

F(t')-n(t')dt> = -J2J2 Xf <y™ / L k{t')L' m {t')-L' k {t')L m {t')dt'. 

k=0 m=0 j 

(41) 

Since LkL' m is a polynomial of order 2 (n — 1) — 1 in t, Gaussian quadrature 
can be used to compute the integral (|4T!) exactly. Finally the integral S/ t+At . 
on the interface at time t + At can be written as 

Sj+l 

S It+AtJ = J F(s') ■ n(s') ds' (42) 

= a t+At (s j+ i) - a t+At (sj), (43) 

where ott+At is the unknown function a representing the interface at time 
t + At. We can solve equation ( 138]) for a t + A t(sj+i): 

&t+At{s j+ i) = M PhPj+l + a t (s j+1 ) - a t (sj) - S Tj + S Tj+1 + a t+At (sj) (44) 

12 



k=0 



= a t (s j+1 ) + S T3+1 - S T0 + M Pk,Pk+i: ( 46 ) 

fc=0 

since we have chosen at(so) = at+i(so) = 0. A kind of similar idea is used 
for surface marker particles in order to correct for area loss during advection, 
cf. 2^1 . In the present method it is, however, not used as a correction but 
as the principle behind the advection step. 
The interpolation points 

(sj, a t+ At(sj)), resp. (s j} Xjyj), j = 0, . . . , N, (47) 

are then interpolated as mentioned in section [3] to give an interpolant of a 
resp. (3 for the interface I(t + At) at time t + At. The area of the drop 
V is conserved, since a(sN = 27r) = V for all time steps. However, this 
conservation property should be understood in a less strict sense, since the 
quadrature errors e a , resp. e& in equations (J2"7j) . resp. fl28l) will lead to the fact 
that the area bounded by the actual interface given by the points computed 
by equations ffTTl) and (|12p can be different from V. In addition, errors can 
lead to the development of self intersections of the interface, cf. figure fl3J). 
Area conservation should rather be understood as the underlying principle 
of the method which is implemented through the function a which, together 
with (3, can be seen as a kind of generating function for the position. 



13 




Figure 2: The fluxing region is the set of points fluxed through the line with end points 
po, resp. pi, from time t to t + At. It is bounded by the trajectories To, resp. n of the 
points po, resp. p\ when tracing them back in time, the line from po to p\ and the image 
of this line when tracing the line back in time. 
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Figure 3: Sketch of the advection step. TOP: Going along the interface I(t) at time t we 
sample a sequence of points pj. MIDDLE: These points are traced forward in time along 
their trajectories Tj to their positions pj at time t + At. At the Gauss Labatto Legendre 
nodes of the interval [t, t + At] , we record the positions of the point pj which are then used 
to find an interpolation approximating the trajectory Tj. BOTTOM: The flux M p . p , +1 
can be interpreted as the signed area of the region bounded by the interface I(t) at time 
t the interface I(t + At) at time t + At and the trajectories Tj and Tj + \. 
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Figure 4: Self intersection of the interface. Errors can lead to self intersection of the 
interface. 
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Figure 5: The error norm E\ measures the total area (marked by red color) included 
between the numerical interface (marked by a blue line) and the exact interface (marked 
by a black line). 



5. Numerical verification 

The numerical verification of the present third order volume of fluid 
method is done by three classical benchmark tests, the reversed single vortex 
test by Rider and Kothe 



tion field test 



19] , Zalesak's slotted disk test [27| and the deforma- 
2l| . However, before going over to the numerical verification, 
we have a glance at the definition of the numerical error. 

5.1. Numerical Error 

The error norm E in the present discussion measures the area difference 
between the numerical interface and the exact interface as depicted in figure 
|5j The error is computed using the position of the interface computed by 
means of equations ( ITT]) , resp. ( TT2j) and not by means of the function a, 
equation ([5]). The integration between the numerical solution and the exact 
interface in order to compute the area difference is done by means of Gaus- 
sian quadrature, where we ensured that the quadrature error is negligible 
compared to the numerical error of the method. The order of convergence O 
between two resolutions n and 2n, is computed using E: 

\n(E(n)/E(2n)) 

~ hT2 • (48) 
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5.2. Numerical verification part 1 



The setup of the reversed single-vortex test of Rider and Kothe [19[ con- 
sists of a circular drop of radius r = 0.15 placed at position (x , y ) = 
(0.5,0.75) in a unit square box. The velocity field u = {—d y ip,d x ip) T is 
obtained by means of the following stream function ip: 

ip(x, y,t) = i cos (jpj sin 2 (nx) sin 2 (ny) , (49) 

where T is the period at which the drop has returned to its initial position. 
Thereby its interface at time T should match the initial circle. The discrep- 
ancy between the numerical interface at time T and the exact circle serves 
as a measure of the numerical error. As initial condition we used 

a = + r (x sin s - y cos s + y ) (50) 

(3 = (r cos s + x ) (r sin s + y ) , (51) 

where s G [0, 2ir\. As mentioned in section |3l we divided [0, 27r] into iV 
equidistant sections. We performed three series of tests for T — 1/2, T — 2 
and T = 8. The results of these tests are shown in figures El d IB1I91 IT01I11II121 
[T3lfT4t [T5lfT6l and tables [U El [3j For all simulations we chose n = 5 for the 
approximation of the trajectories for the advection step, cf. section HI De- 
pending on the order P of the B-spline interpolation we chose a different 
value for the time steps At in order to make the error contribution due to 
the advection step subdominant compared to the error contribution of the 
interface Reconstruction. The advection step itself can handle quite large 
time steps without displaying any sign of instability. All simulations were 
performed using B-splines of order P = 3,5, 7. For P = 3 the time step was 
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chosen At = 1/32, for P = 5, At = 1/128, and for P = 7, At = 1/256. 
The time step was kept fixed at these values even when going over to finer 
resolutions. In figure [6] the resulting position of the interface is shown for 
t = T/2 and t = T in the case T = 1/2, P = 3 and a resolution of N = 10. 
The interface is well resolved and does not display any visual disturbances 
such as bumps or oscillations at maximum deformation and when it has re- 
turned to its initial position at t = T. The same observation can be made for 
the case T = 2, cf. figure El Concerning the convergence of the method for 
these two cases, e.g. T = 1/2 and T = 2, cf. figures El resp. [9] and tables HJ 
resp. [21 we observe that the order of convergence corresponds approximately 
to the theoretical value of P + 1, cf. equation ([2l]) . For P = 7 and fine 
resolutions the error does not further decrease because of the round off limit. 
In addition, for the case T = 2 when going from iV = 20 to iV = 40, we 
observe a sudden jump in the convergence, cf. figure [9] and table [21 for P = 5 
and P = 7. This phenomenon of accelerated convergence is even more pro- 
nounced in the case T = 8, cf. figure [121 and table [31 An explanation for this 
sudden increase in convergence might lie in a underresolution of the problem 
for coarse resolutions N, meaning that when increasing P, keeping N fixed 
at a small value no important gain in accuracy is observed. However, if the 
resolution is finer, i.e. large values of N, increasing P will almost lead to 
spectral convergence, i.e. faster than algebraic. This points to the eventual- 
ity of having insufficient sampling of the signal, i.e. an aliasing error. This is 
different to classical spectral methods such as methods based on Chebyshev 
or Legendre polynomials for which increasing the approximation order in- 
troduces an increase of spatial resolution by increasing the number of Gauss 
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points. Periodic B-splines on the contrary offer the possibility of increasing P 
and N independently with the consequence of having eventually a persisting 
aliasing error when only increasing P. Even worse a underresolved simula- 
tion can lead to an entirely wrong solution by the present method, as can be 
seen in figure [TUl Here T = 8 and the resolution N was fixed at 20 (P = 3). 
For t = T/2 still a few structures of the correct solution are recognizable, 
however for t = T the picture has entirely deteriorated. If underresolved, 
the numerical solution can display self intersections. This indicates that the 
present method is less robust to underresolution. However, for well resolved 
cases, which for the case T = 8 start at only 40 knots, cf. figure HH the 
present numerical scheme produces a very accurate result. 
Concerning area conservation, we observe from figure [12] that the value a(27r) 
is, apart from round off contributions, equal to the initial area of the drop. 
However, the quadrature errors e a , resp. e&, equations (|27|) . resp. ( 128|) . can 
become rather important during simulation, as for instance for the underre- 
solved case T = 8, iV = 20 and P = 3, cf. figure HU indicating a possible 
discrepancy between the actual area of the drop and a(27t). For well resolved 
cases the quadrature errors are smaller, cf. figure HU but seem to increase 
with increasing deformation of the drop. 

As a last numerical experiment we investigated the accuracy of the advec- 
tion scheme derived in section [H By fixing the number of Gauss Labatto 
Legendre nodes to n = 5, the interpolating polynomial has order 4 for which 
reason we expect the advection scheme to converge with 5 th order accuracy 
with respect to the time step At. In order to observe the error contribution 
by the advection scheme, we chose a B-spline of order P = 7 and a spa- 



20 




Figure 6: Result of the Rider-Kothe single vortex benchmark test for T = 1/2. The 
resolution was N = 10 and P = 3. LEFT: Position of the interface at t = T/2. RIGHT: 
Position of the interface at t = T. 

tial resolution of N — 160, such that the error contribution by the interface 
representation is subdominant. Decreasing the time step At leads indeed to 
a fifth order convergence of the numerical error up to the point at which 
the error contribution by the interface representation becomes dominant, cf. 
figure [16j 
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P = 3 
slope -4 

P = 5 
slope -6 

P = 7 
slope -8 



Figure 7: Error decrease for the Rider-Kothe single vortex benchmark test for T = 1/2, 
for increasing resolution N using B-splines of different order P. 




Figure 8: Result of the Rider-Kothe single vortex benchmark test for T = 2. The resolution 
was N = 10 and P = 3. LEFT: Position of the interface at t = T/2. RIGHT: Position of 
the interface at t = T. 
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Figure 9: Error decrease for the Rider-Kothe single vortex benchmark test for T = 2, for 
increasing resolution N using B-splines of different order P. 




Figure 10: Result of the Rider-Kothe single vortex benchmark test for T = 8. The 
resolution was N — 20 and P = 3. LEFT: Position of the interface at t = T/2. RIGHT: 
Position of the interface at t = T. In this case the resolution chosen was too coarse leading 
to a break down of the present method. 
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Figure 11: Result of the Rider-Kothe single vortex benchmark test for T — 8. The 
resolution was N = 40 and P = 3. LEFT: Position of the interface at t = T/2. RIGHT: 
Position of the interface at t = T. 
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Figure 12: Error decrease for the Rider-Kothe single vortex benchmark test for T 
increasing resolution N using B-splines of different order P. 
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.R = 3, N -20 
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Figure 13: The difference of the exact area V of the drop to the area of the drop given by 
a(2ir) during the Rider-Kothe single vortex benchmark test for T = 8. 
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Figure 14: The quadrature errors e a , resp. e&, defined in equations (|27|) . resp. ([28]) during 
the Rider-Kothe single vortex benchmark test for T = 8. The corresponding position of 
the interface is displayed in figure \W\ at some instances in time. For this low resolution 
(N = 20, P = 3) the method breaks down. This is also indicated by large quadrature 
errors meaning that the actual position of the interface, computed by means of equations 
(fTTjl . resp. (fl2]l. includes a large deviation from the exact position given by a and (3. 



25 



e_, P = 3, N = 40 
E b , P = 3, N = 40 
e,, P = 3, N = 80 
E b , P . 3, N = 80 



150 

number of iterations 



Figure 15: The quadrature errors e Q , resp. eb, defined in equations (|27[) . resp. (|28|) 
during the Rider-Kothe single vortex benchmark test for T = 8. This time the resolution 
is sufficient (N = 40,80, P = 3) which is also indicated by a smaller quadrature error 
compared to the one in figure [T4l For maximum deformation of the drop we are confronted 
with a maximum quadrature error. 
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Figure 16: Error increase for the Rider-Kothe single vortex benchmark test for T = 8, for 
increasing time step At using B-splines of order P = 7 and a resolution of N — 160. 
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Table 1: Results for the Rider-Kothe single vortex benchmark test for T 



p 


N 


Error 





3 
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2.96591 x 10- 


-4 






20 


3.54562 x 10" 
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2.51014 x 10- 
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Table 2: Results for the Rider-Kothe single vortex benchmark test for T 
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p 


N 
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2.34006 x 10- 


-5 






80 


4.21891 x 10- 


-8 
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3.82007 x 10" 


-10 


6.79 
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5.19595 x 10- 


-12 


6.20 
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40 


1.16624 x 10- 


-5 






80 


7.93665 x 10- 


-9 


10.52 
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1.92373 x 10- 


-12 


12.01 




320 


3.63856 x 10- 


-13 


2.40 



Table 3: Results for the Ridcr-Kothe single vortex benchmark test for T = 8 
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5. 3. Numerical verification part 2 

The slotted disk rotation test of Zalesak |27] uses a. solid body rotation 
to advect a slotted disk. The stream function ip is given by 

^{x,y) = ^{(x-x ) 2 + (y-y ) 2 }, (52) 

where u is chosen in such a way as to allow a complete rotation of the drop 
in 2524 iterations. The computational box is four by four and the drop has 
a radius of r = 1/2. It is situated at (x ,y ) = (2,2.75). The coordinates of 
the four corners of the slot are given by: 

(xa,Va) = (x + ^,y ~ ]frj- (x a - x ) 2 ) (53) 

5 

(x b ,y b ) = (x a ,yo + -^) (54) 
3 

(xc,Vc) = (x -—,y b ) (55) 

(,Xd,yd) = (x c ,ya)- (56) 

As an initial condition we chose a description of the initial interface by means 
of four functions a±, a 2 , a 3 , a 4 , defined the following way: 

Qfi(s) = - (r^s - r (x (cos (s + s a ) - cos (s a )) + y (sin(s + s a ) -sin(s a 
a 2 (s) = ^x c s, (58) 

«s(s) = ~\vb s , ( 59 ) 
a 4 (s) = -x a s, (60) 

(61) 
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where s a = arcsin ' Xa xo 



( Xa ro X ° ) anc ^ Sb = arcs i n ) • The function a is com- 



posed by means of these four functions. 

an(s) 

a 2 {s - 2tt + s b + s a ) 

+«l(27r - S 6 - S a ) 

a 3 (s - 2tt + s 6 + s a - y 6 + y a ) 

+a2(Vb - Da) 

+ai(2n - s b - s a ) 
a 4 (2vr - s b - s a + y b - y a + x a - x c - s) 
+a 3 (a;a - £ c ) 

+ttl(27T - S 6 - S a ) 



a(s) 



< S < 27T — s b — s a 

2lT - S b - S a 

< S < 27T - S b - S a + 2/6 - 2/a 

2tt - s 6 -s a + y b -y a 

< s 

< 2tt - Sb - s a + y b - y a + x a 

2n- s b - s a + y b -y a + x a - x, 

< s 

< 2ir - Sb - s a + 2(yb - y a ) + 

(62) 

- Sb - s a + 2(y 6 - y a ) + 



The parameter s takes values in the interval [0, 2n 
x a — x c] this time. We discretized this interval in such a way that the corner 
positions of the slot are at knots of the discretization. The function a is then 
interpolated at these knots. The function (3 is handled likewise, with j3 given 
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by: 



[r sm(s + s a ) + xo) 
(-r cos(s + s ) +y ) 



< s < 2n — s b — s a 



Xd{y a + s - 2n + s b + s a ) 
Vb{xd + s - 2tt + s b + 

s a ~ Vb + Va + Xd) 

Xa{Vb - S + 2-K - S b - 
Sa + Vb-Va + Xa- X c ) 



2-k - s h - S a 

< s < 2tx - s b - s a + y b - y a 
2tt - s b - s a + y b - y a 

< s 

< 2ir - s b - s a + y b - y a + x a - x c 
27r - s b - s a + y b - y a + x a - x c 

< s 

< 2tt - s b - s a + 2(y b - y a ) + x a - x c . 

(63) 

The initial condition is only C° because of the kinks at the corners of the slot. 
Although using high order periodic B-splines, we expect the convergence rate 
therefore to be of only second order at most. In addition, discontinuities can 
give rise to the development of spurious oscillations of the interpolant at these 
discontinuities, the so called Gibbs phenomenon [5]. The higher the order 
of the periodic B-splines the further these spurious oscillations spread along 
the interface, as can be seen when comparing figures [T7| and [TS], where we 
compare the interface position at time t — and t = T, after one rotation, for 
different resolutions. Intermediate steps are shown for iV = 160 and P = 3 
in figure [HE indicating that the main contribution to the overall error does 
indeed not come from the advection but from the interpolation of a function 
with discontinuities in its first derivative. The order of convergence is reduced 
to second order no matter the order of the B-spline interpolation, cf. figure 
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1.58615 x 10" 


-2 
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4.25359 x 10" 


-3 


1.90 
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1.08273 x 10- 


-3 


1.94 


5 


28 


3.76212 x 10° 








56 


2.40403 x 10" 


-1 


3.97 




112 


4.16235 x 10- 


-2 


2.53 




224 


9.12356 x 10- 


-3 


2.19 




448 


2.12014 x 10- 


-3 


2.11 



Table 4: Results for Zalesak's slotted disk rotation test for P — 3 and P = 5. 

[20| resp. table HI The absolute error is even larger for higher order B- 
spline interpolation, due to the larger spurious oscillations. We remark that 
choosing an order P = 1 for the B-spline interpolation gives us a PLIC like 
description of the interface, which for the present benchmark test produces 
more accurate results since it only requires C° continuity of the function to 
be interpolated. This leads to the result that for P = 1 the slot stays sharp 
during the entire simulation as can be seen in figure [2TJ 
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jV = 28,t = N = 28,t = T 




N = 448, t = N = 448, t = T 



Figure 17: Result of Zalesak's slotted disk rotation test for different resolutions N. The 
order of B-splinc interpolation is P = 3. Shown are the graphs at the initial position, 
t = 0, and the position after a full rotation, t = T. 
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N = 28, t = 



N = 28, t = T 




N = 448, t = N = 448, t = T 



Figure 18: Result of Zalesak's slotted disk rotation test for different resolutions N. The 
order of B-splinc interpolation is P = 5. Shown are the graphs at the initial position, 
t = 0, and the position after a full rotation, t = T. 
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t = 



t = T/4 





Figure 19: Result of Zalesak's slotted disk rotation test for N — 224 and P = 3. The 
graphs are shown at different points in time. 
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Figure 20: Error decrease for Zalesak's slotted disk rotation test. Due to the discontinuity 
in the first derivative of the solution, the present method exhibits only second order accu- 
racy with respect to the resolution N. In addition, the fifth order B-splinc interpolation 
is less accurate than the third order one, since the spurious oscillations are larger in the 
former case. 
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t = 



t = T/4 




t = T 



Figure 21: Result of Zalesak's slotted disk rotation test for N = 28 and P = 1. The 
graphs are shown at different points in time. 
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5.4- Numerical verification part 3 
The deformation field test 



21] uses the following stream function: 



ij>{x, y, t) = ^ cos ^jpj sin ^nn (% + 7^jj cos (^ n7r \ V + J > ( 64 ) 
where n is the number of vortexes in the computation domain and chosen 



to be n = 4 to match the geometry used in [28|, |24j|, as was the period with 
T = 2. The computational domain is a square box of side length one. A drop 
of radius r = 0.15 has its center at (xq, yo) = (0.5, 0.5) at time t — 0. Since 
the flow is reversed after T/2 the drop returns to its initial position at t — T 
assuming its initial shape. As initial condition we used the same functions 
a, resp. (3, as for the Rider-Kothe single vortex test, equations fl50H5ip . We 
performed a series of tests with P = 1, P = 3 and P = 5. The time steps 
were chosen At = 3.906 x 10" 3 for P = 1 and P = 3 and At = 9.76 x 10~ 4 for 
P = 4. This was, as before, done in order to make the error contribution by 
the advection step subdominant. This benchmark test is relatively difficult 
since the interface develops regions with very small local radii of curvature. 
In addition, at some parts the drop becomes very thin. The results of the 
present method applied on this benchmark test can be seen in figures [22] for 
P = 1, [23] for P = 3, and [241 for P = 5. In these figures we depicted the 
graph of the interface at maximum deformation t = T/2 and after the drop 
has returned to its initial position at t = T for different resolutions. For low 
resolutions the graph at maximum deformation (t = T/2) does only capture 
the coarse features of the solution for all three orders P = 1, P = 3 and 
P = 5. In addition, at regions were the resolution is low but the curvature 
high, the numerical solution seems to develop a kind of Gibbs phenomenon 
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for P = 3 and P = 5, as for Zalesak's slotted disk rotation test. These 
oscillations become smaller as the resolution increases. When the drop has 
returned to its initial position we observe that the interface develops spikes 
at those regions where the resolution was low at maximum deformation. This 
is due to the fact that as the error is larger in these regions a point might 
fall into the wrong vortex and be traced to a different location. These spikes 
become smaller with increasing resolution N and increasing order of the 
interpolating B-spline P. For finer resolutions, N = 1250, the final position 
is extremely close to the exact solution. However, at maximum deformation 
some wiggles can be observed for both P = 3 and P = 5, indicating that 
measuring the error at maximum deformation might give a better estimate 
of the accuracy of the present method than measuring it at the final position. 
For iV = 5000 and P = 5 the wiggles have disappeared, as can be observed 
from figure [251 Concerning the order of convergence, cf. figure [26] and table 
the present method seems to converge at a lower speed both for P = 3 
and for P = 5. This might have its origin in the spurious oscillations which 
might prevent the method from converging at the right rate. Grid adaptation 
or remeshing redistributing the points to regions were needed, as mentioned 
in section [31 might be advantageous for this benchmark test. Nevertheless, 
an order of convergence between three and four with respect to the spatial 
resolution for this benchmark test is quite acceptable. 



39 



» 


[T 


p 




N = 


156, t = T/2 


„ 


i 


1 

23 

■ ■ . -fJ 




N = 


312, i = T/2 


„ 








bp 




N = 


625, t = T/2 




b 

V 


--™«*^| 
> 





AT = 1250, t = T/2 




AT = 156, t = T 




N = 625, t = T 




N = 1250, t = T 



Figure 22: Result of the deformation field test for different resolutions N. The order 
of B-spline interpolation is P = 1. The graphs are shown at the maximal deformation, 
t = T/2 and after returning to the initial position t = T. The black dashed line is the 
exact solution for t = T. 
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N = 156, t = T/2 



N = 156, t = T 






Figure 23: Result of the deformation field test for different resolutions N. The order 
of B-spline interpolation is P = 3. The graphs are shown at the maximal deformation, 
t = T/2 and after returning to the initial position t = T. The black dashed line is the 
exact solution for t = T. 
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N = 156, t = T/2 N = 156, t = T 





Figure 24: Result of the deformation field test for different resolutions N. The order 
of B-spline interpolation is P = 5. The graphs are shown at the maximal deformation, 
t = T/2 and after returning to the initial position t = T. The black dashed line is the 
exact solution for t = T. 
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Figure 25: Results of the present method for the deformation field test at different points 
in time, here N = 5000 and P = 5. 
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5.30438 x 10" 
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3.33 
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1.11787 x 10- 


-3 
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1.14428 x 10- 


-4 


3.29 




625 


1.12848 x 10- 


-5 


3.34 
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1.40122 x 10- 


-6 


3.01 
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1.07634 x 10- 


-7 


3.70 
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1.14465 x 10- 


-8 


3.23 



Table 5: Results for deformation field test for P = 1,P = 3 and P = 5. 
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Figure 26: Error decrease for the deformation field test. The poorer convergence might 
be due to the appearance of spurious oscillations at regions of low resolution. 
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6. Conclusions 

In the present discussion we derived an alternative formulation for the 
interface representation for the volume of fluid method. The interface is 
represented in a periodic fashion by two functions a and j3, from which the 
position of the interface can be calculated. These two functions are approxi- 
mated by periodic B-spline interpolation which allows a systematic extension 
to higher order accuracy with respect to the grid spacing. The advection 
scheme has been simplified and extended to higher order accuracy with re- 
spect to the time step. Numerical verification indicates that the present 
scheme has indeed the order of convergence predicted by the theory. This 
allows for very accurate simulations with a limited number of knots. How- 
ever, if the sampling rate is too small, the present scheme can break down, 
providing a numerical solution far away from the exact one. In addition, 
for discontinuities in the first derivatives at a point on the interface, or at 
regions of the interface with poor resolution and high curvature, the present 
method can display a kind of Gibbs phenomenon. Taking a lower order B- 
spline interpolation P = 1 can in this case provide more appealing results. 
A remeshing or adaptive grid approach could increase the efficiency of the 
algorithm. In addition, such an approach might furnish a way to simulate 
topological changes such as coalescence or drop break up. The extension to 
three dimensions is also left for future research. 
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